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ABSTRACT 

This  paper  derives  the  transport  equations  for  rarefied  gases  from  the  Bhatnagar-Gross-Krook  (BGK)  model 
kinetic  equation  using  Hermite  polynomial  representation  of  the  velocity  distribution  function.  We  apply  the 
Champmen-Enskog  method  to  Grad’s  thirteen  moment  equations  to  derive  a  closure  of  Grad’s  13  moment 
equations,  extending  them  to  third  order  of  the  Knudsen  number.  The  velocity  distribution  function  for  the 
resulting  13  regularized  moment  equations  is  presented. 


I.  INTRODUCTION 

One  of  the  hardest  problems  in  computational  fluid  dynamics  is  the  modeling  of  medium  rarefied  gases 
with  the  Knudsen  number  in  the  range  0.005  -  1.  In  this  case,  the  gas  is  rarefied  to  such  a  degree  that  using  the 
Navier-Stokes-Fourier  equations  is  questionable,  but  it  is  not  sufficiently  rarefied  that  using  Direct  Simulation 
Monte  Carlo  (DSMC)  methods  is  effective.  It  should  be  stressed  that  in  recent  years  computational  fluid 
dynamics  for  Knudsen  numbers  in  this  range  has  become  more  and  more  important  for  many  practical 
applications,  ranging  from  the  modeling  of  reentry  of  space  vehicles  into  the  atmosphere  to  the  modeling  of 
microscale  flows  and  heat  transfer  in  microchannels.  One  of  the  ways  to  cover  this  range  of  Knudsen  numbers 
is  to  use  thirteen  (or  more)  moment  equations  instead  of  the  Navier-Stokes-Fourier  5  equations.  In  1949  Grad 
derived  13  moment  equations  corresponding  to  the  second  order  of  the  Knudsen  number  [1,  2].  It  is  worth 
noting  that  the  Navier-Stokes-Fourier  equations  are  to  first  order  of  the  Knudsen  number.  Unfortunately,  Grad’s 
moment  equations  sometimes  produce  unphysical  solutions;  for  example,  they  fail  to  describe  smooth  shock 
structures  for  Mach  numbers  above  a  critical  value  [3].  In  2004  Struchtrup  [4]  regularized  Grad’s  equations, 
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extending  them  to  third  order  of  the  Knudsen  number.  This  was  a  very  important  step  in  rarefied  gas  dynamics. 
This  step’s  significance  can  be  compared  to  the  extension  of  Euler’s  gas  dynamics  equations  to  the  Fourier- 
Navier-Stokes  gas  dynamics  equations.  The  author  has  developed  a  new  closure  method  that  is  principally 
different  from  the  well-known  Chapmen-Enskog  method  [5,  6]  (that  was  used  to  derive  a  closure  of  Euler’s  gas 
dynamics  equations),  in  which  he  has  not  used  the  Hermite  polynomial  representation  of  the  velocity 
distribution  function.  It  should  be  noted  that  the  Struchtrup  closure  method  is  very  complicated,  and  this  is 
probably  one  of  the  reasons  that  his  method  and  equations  are  difficult  to  comprehend. 

In  the  present  paper  we  suggest  a  new  closure  for  Grad’s  thirteen  equations  by  using  a  Hermite  polynomial 
approximation  for  the  monatomic  gas  velocity  distribution  function,  and  applying  the  Chapman-Enskog 
regularization  method  to  Grad’s  velocity  distribution  function  that  corresponds  to  his  13  moment  equation  [5, 
6].  In  our  paper,  the  collision  term  is  assumed  to  be  in  the  BGK  form.  The  integral  representation  for  the  13 
moments  of  the  Boltzmann  equation  and  the  Hermite  polynomial  approximation  of  the  velocity  distribution 
function  are  obtained  in  Sections  II  and  III  respectively.  Grad’s  regularized  13  moment  equations  are  derived  in 
Section  IV,  and  conclusions  are  presented  in  Section  V. 


II.  GENERAL  EQUATION  FOR  13  MOMENTS 


The  phase  density  of  a  monatomic  ideal  gas  is  described  by  the  Boltzmann  equation. 


d{n-f)  ^  d(n-f) 
dt  *  dxi 


=  St(n-f), 


where 


V 


(1) 

(2) 


Here  n  is  the  number  density  of  gas  molecules,  /  is  the  velocity  distribution  function,  Vi=(V„Vy,V,)  is  the 

particle  velocity,  Xj  =  (x,y,z)  are  the  coordinates  of  a  particle,  j  d^V  means  the  integration  over  the  entire  velocity 

V 

space,  and  St(n  ■  f)  is  the  collision  term  that  accounts  for  the  change  in  the  velocity  distribution  function  due  to 
collisions.  Here  we  assume  elastic  collisions. 

Let  us  introduce  p  as  the  mass  density  of  gas  molecules. 


P 


V 


(3) 
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Uj  =  (u^,Uy,u^)  as  the  flow  velocity  of  gas  molecules, 

ui^lf-Vrd^V,  (4) 

V 

Vj  as  the  thermal  velocity, 

qi  =  {qx  ’  the  heat  flux. 


qi=^-lf{V)-{Vi-Ui)Wx-UxY+iVy-UyY+{V^-u,Y]-d^V, 


\2i  j3i 


(jy  =  {<j^,<j^y,<j^2,ayy,<jy^)  as  thc  compoucnts  of  the  stress  tensor. 


^ii  =P-\f 


7/2  ^ 

{Vi-Ui)-{Vj-Uj)-dy-^ 


■d^V, 


(6) 


(7) 


were  m  is  the  mass  of  a  particle.  As  one  can  see,  the  total  number  of  moments  of  the  particle  distribution  function 
introduced  here  is  13.  Let  us  obtain  general  equations  for  these  13  moments  using  the  Boltzmann  equation,  Eq. 

(1).  First  we  consider  equations  for  p  ,  ,  Uy,  ,  and  .  Since  the  number  of  colliding  particles,  their  total 

momentum,  and  their  total  energy  are  conserved  in  collisions,  it  follows  that 

fi-fSt(f)-d¥^0,  (8) 

V 

m-n-\VpSt(f)-d^V  =  0,  (9) 

V 


—  ■n-\V^-St(f)-d^V  =  0.  (10) 

2  f 

We  obtain  from  Eq.  (1)  the  following  moment  equations  that  correspond  to  mass,  momentum  and  energy 
conservation  laws  respectively. 


dt 


m-n-  \  f-d^V 
V  J 


dX; 


m-n  -  \  f  -  Vi  -  dy 
\  V  J 


-0, 


(11) 
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—  m-n- 1  f -Vj^ -d^V  +^—  m-n- \  f d^V  =0, 

dt\  y  )  dXi\  y  ) 


dt\  2 


dXi  \  2  f 


where  =  V^  +Vy  +  V^  ,  with  indexes  i,k  =  x,y,z  .  After  tedious  algebra  Eq.  (1 1)  -  (13)  can  be  presented  in  the 
following  form  [7] 


dp  d  ,  . 


duj  d  { p '  Vf  ^  d(^ji 
— -  +  —  - — —  + — ^  =  0, 


p-—^  +  (p-u,)-—^  +  —  ^ 

dt  ^  dx  j  dx:  2  dx ; 

j  I  \  y  j 


3  dVf  3  0  (  Fj  1  1  2  du:  dq;  du; 

-■  p - —  +  -■  p-Uj - —  +-■  p-Vt  ■  — -  +  —  +  (T,7 - -  =  0  . 

4  dt  2  dxj  I  2  J  2  dxj  dxj  ^  dx j 


We  derive  the  general  moment  equations  for  and  a^y  and  in  Appendix  A.  They  can  be  written  as: 


ddi  I  5  P'Vt]  duj  duu  dr  x  du;  (  5  p-Vy]  dw  /  ^  du.- 
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where  =  cjji ,  {V -uf  =  (F^  -u^f  +{Vy  -Uyf  +(F^  -u^f ,  and 


^zz  =  p-\  f -^-^=-^yy-^zz-  (20) 

V  ^ 

Changing  the  order  of  indexes  {x,y,z)  to  {x,z,y)  in  Eq.  (18)  we  obtain  an  equation  for  , 


da^ 


dUy 


5m, 


a,  ■  +  - +  (^kx'^ - +  ^ — 

Ot  UX  UX  ^  UX 
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=  p-\{V,-u,)-{V,-u,)-St{f)-d^V. 


and,  proceeding  similarly  with  (z,y,x)  we  obtain  an  equation  for  ^ 
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ZF 
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+  (^ky  •  ^ ^ ^ — 
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Changing  the  order  of  indexes  (x,y,z)  to  (y,x,z)  in  Eq.  (19)  we  obtain  an  equation  for  o 


yy 


1  t^2  2  dq;  2  du;  d  (  \  ^  ^^2  ^^y 
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It  is  worth  noting  that  there  are  no  collision  terms  in  the  moment  equations  for  p  ,  ,  and  F^ ,  Eqs.  (14)  - 

(16),  while  the  equations  for  heat  flux  and  stress  tensor,  Eqs.  (17)  -  (19)  and  (21)  -  (23)  do  include  collision  terms. 
Since  collision  terms  can  produce  new  moments,  the  set  of  these  13  moment  equations  in  general  is  not  self- 
contained.  However,  in  the  case  of  Maxwell  molecules  and  the  BGK  approximation  of  the  collision  term, 

St{f)  =  ^^^^,  (24) 

T 


5 


where  r  is  a  collision  time  depending  on  time  and  coordinates  and  is  the  Maxwellian  velocity  distribution 
function, 


fn  - 


r  ^  a3/2 


yn-Vj  j 


•exp 


^  {V,-U,f  +  {Vy-Uyf  + 

Vt 


(25) 


the  collision  terms  do  not  produce  any  new  moments.  In  other  words,  these  approximations  of  the  collision  term 
do  not  mix  the  moments.  This  is  a  key  point  for  any  theory  of  moment  approximation  of  the  Boltzmann  equation. 
In  this  paper  we  will  use  the  BGK  approximation  of  the  collision  term;  the  case  of  Maxwellian  molecules  can  be 
described  in  a  similar  way. 


III.  HERMITE  POLYNOMIAL  APPROXIMATION  OF  FUNCTION  DISTRIBUTION 

Following  Grad  [1,  2]  we  assume  here  a  Hermite  polynomial  approximation  of  the  velocity  distribution 
function  where  Hermite  polynomials  are  described  as  follows: 

Hn{x)  =  (-1)^  •  exp(j^ )•  [exp(-  )]  ,  (26) 

dx 

The  velocity  distribution  function  can  be  described  as  a  combination  of  three-dimensional  Hermite  polynomials 
that  correspond  to  x ,  y  ,  z  directions  of  the  velocity.  For  our  purposes  we  need  only  the  following  set  of  Hermite 
polynomials: 

«0=1.  ■Wlito)  =  2'Zi.  W2,te)  =  4-^?-2,  (28) 

ff3i(x)  =  8'2'?-12-2'i,  ff4i(2')  =  16-x;'-48-x?  +  12,  (29) 

where  i  =  x,y,z  and 
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(30) 


Vj-Uj 

Vj 


We  represent  the  velocity  distribution  via  the  29  Hermite  polynomials 


29  .  ,  . 

fa  /m  '  ( Xx  ’  Xy  5  Xz  )  ’ 

k=\  \  X  ; 


(31) 


where  4^  is  a  Maxwellian  function,  Eq.  (25),  coefficients  A  depend  on  coordinates  and  time,  and  Hermite 
polynomials  H  are 
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It  should  be  stressed  that  all  H  polynomials  are  orthogonal,  i.e., 


2  .2  2 


J  H^-Hj-exp[-Zx-Zy-Zz]-dz  =  0,  where  i^j 


lx  yiy  yiz 


(48) 


Since  the  velocity  distribution  function  has  to  satisfy  the  conditions  given  by  Eqs.  (2),  (4),  and  (5),  we  obtain 
that 


2^5  +2Vg  +  Ajo  -  0  . 


=  1 ,  A2  ^  A3  ^  A4  ^  0  and 


(49) 


Thus,  the  particle  distribution  function  n-f^  has  29  variables,  p ,  u^,  Uy,  ,  Vj  ,  A5 ,  Ag,  A7,  Ag ,  A9 ,  Ajj- 
A29  .  Substituting  the  velocity  distribution  function  ,  Eq.  (31),  for  /  in  Eqs.  (6)  and  (7),  we  obtain  relationships 
between  ,  ay  and  A  : 
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It  is  worth  noting  that  the  truncated  velocity  distribution  function  that  consists  of  the  first  nine  nonzero  Hermite 
polynomials  has  the  form  of  the  Champmen-Enskog  and  Grad’s  velocity  distribution  functions  [1,  2,  7].  In  the  next 
sections  it  will  be  shown  why  we  have  selected  this  representation  of  velocity  distribution  function. 


IV.  A  CLOSURE  OF  GRAD’S  13  MOMENT  EQUATIONS. 

Let  us  rewrite  the  equation  for  the  heat  flux,  Eq.  (17),  and  stress  tensor,  Eqs.  (18)  -  (19),  (21)  -  (23)  for  the 
case  of  the  BGK  collision  term,  Eq.  (24): 
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where  indexes  i,k  =  x,y,z  .  To  complete  this  system  of  equations,  a  velocity  distribution  function  has  to  be  chosen. 
Grad  [1,2]  has  suggested  the  following  velocity  distribution  function: 

foRAoi^)  =  fu 

and  obtained  his  set  of  equations  that  later  has  been  called  “Grad’s  13  moment  equations.”  As  one  can  see  Grad’s 
velocity  distribution  function  can  be  rewritten  in  terms  of  Hermite  polynomials  as 
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where  we  have  taken  into  account  that  oy  =o-ji  and  used  Eq.  (20).  Grad’s  equations  [1,2]  for  and 

can  be  presented  in  the  following  form: 
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d  /  \  du^ 


^  P-Vt^ 
V  2  , 


0M„  SUy 

— ^  +  <^Ax  • — -  + 
dy  dxk 


^ p-Vj^  Su 


V  2  , 


V  ^  2  ^  2  dq^  _  <^xy 


dx  5  dx  5  dy 


Changing  the  order  of  indexes  (x,y,z)  to  (y,x,z)  and  (z,y,x)  in  Eq.  (59)  we  obtain  the  corresponding  equations 
for  qy  and  ;  likewise,  using  (y,x,z)  in  Eq.  (60)  we  obtain  an  equation  for  <jyy-,  and  finally,  using  (x,z,y) 
and(y,z,x)  in  Eq.  (61)  we  obtain  equations  for  and  respectively.  Unfortunately,  as  it  has  been  mentioned 

in  Introduction,  Grad’s  moment  equations  sometimes  produce  unphysical  solutions  and,  therefore  need  to  be 
regularized. 

Let  us  obtain  a  set  of  the  13  regularized  Grad’s  moment  equations  using  the  Hermite  polynomial 
approximation  of  the  velocity  distribution  function  and  the  Chapmen-Enskog  closure  method,  while  assuming  that 
the  collision  term  is  in  BGK  form,  Eq.  (24). 

First  let  us  represent  the  polynomials  in  Eq.  (51)  with  index  i  =  x  and  in  Eqs.  (52),  Eq.  (55)  in  Hyrmite 

forms: 
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2  1)  ..  ^3x  ,  ^Ix 
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Zx  2  1  8 


(65) 


(66) 


(67) 


(68) 


(69) 


(70) 


(71) 


(72) 


where  variables  j  are  given  by  Eq.  (29).  As  one  can  see,  the  polynomials  presented  in:  Eqs.  (62)  -  (68)  are 
included  in  ;  Eqs.  (66),  (68)  and  (69)  in  cr^y ;  and  Eq.  (70)  -  (72)  in  a^x  ■  The  Hermite  representations  of 
polynomials  in  Eq.  5 1  with  index  i  =  y,z ,  and  in  Eqs.  (53),  (54),  (56)  for  qy,  q^,  ^  <^>'z  ^  ^>^0  <7yy  respectively 

is  obtained  from  Eqs.  (62)  -  (72)  by  the  proper  rotation  of  the  indexes.  Thus,  the  complete  list  of  Hermite 
polynomials  that  represents  polynomials  in  Eqs.  (5 1)  -  (56)  is 

//4x,  H^y,  (73) 
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(74) 


^2x  '  ^2y  ’  ^2x  '  ^2z  ’ 


ti2y-H2., 


^3x  ’  ^ly  1  ^3x  '  ^Iz  1  ^3y  '  ^Ix  ’ 


Hlx  ■  H,y  ■  //2z  ,  iixx  ■  ^Iz  ■  ^2,. ,  iixz  '  ^Xy  '  ,  (76) 

H^x-Hxy-Hy,,  ill) 

7/3,,  7/33,,  773,,  (78) 

^2x  ■  ^Xy  5  ^2x  ■  ^Xz  ■'  ^2y  '  ^Xx  ■'  ^2y  '  ^Xz  5  672^  ■  77^^  ,  772^  ■  Hyy  ,  (79) 

H2x^  f^2y,  672,,  (80) 

771.. 7713.,  77i,.77i,,  77i3,.77i,,  (81) 

661.,  77i3,,  77i,,  (82) 

77o.  (83) 

Now,  following  the  Chapman-Enskog  method  [5,  6],  let  us  write  the  velocity  distribution  function  as 
/"  ^  fcRAD  +  ^  ■  /m  ■  /l  •  (84) 


where  has  a  Hermite  polynomial  form.  Since  the  velocity  distribution  function  t-  fy^-  does  not  have  to 
contribute  into  the  previously  obtained  13  moments  p,  w, ,  and  Vj  ,  ,  and  ap ,  it  follows  that  the  Hermite 

polynomials  included  in  Grad’s  velocity  distribution  function,  Eq.  (58),  have  to  be  excluded  from  fj .  But  as  the 
velocity  distribution  function  r- /m'  fx  6as  to  contribute  into  the  integrals  in  Eqs.  (51)  -  (56),  we  obtain  that  /J , 
as  follows  from  Eqs.  (62)  -  (72),  has  to  be  a  combination  only  of  the  Hermite  polynomials  presented  in  Eqs.  (73)  - 
(77).  Subsequently,  we  obtain  that 

29 

/i=ZA-77,  ,  (85) 

i=XA 
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where  functions  i/,  are  shown  in  Eqs.  (43)  -  (47).  Thus,  we  have  shown  that  the  chosen  set  of  Hermite 
polynomials,  Eq.  (32)  -  (47)  has  good  physical  sense  for  representing  the  velocity  distribution  function. 

Eet  us  introduce  the  following  set  of  M-moments: 


V 

M2aj-V^-p-\f-H2rH2j-d^V, 

V 

V 

Muijlk  =  VT-p-if-Hu-  Hlj  ■  Hit  ■  V  , 

V 

Muijlt  =  Vhp- Lf-Hu-Hij -Hit -d^V  , 

V 


i  =  ix,y,z),  (86) 

ij  =  {xy,?z,yz) ,  (87) 

ij  =  {xy,xz,yx,yz,zx,zy) ,  (88) 

ijk  =  {xyz,yzx,zxy),  (89) 

ijk  =  (xyz),  (90) 


which  correspond  to  Hermite  polynomials  in  Eqs.  (73)  -  (77).  Substituting  them  into  Eq.  (51)  with  index  i  =  x  and 
into  Eqs.  (52)  and  (55)  we  obtain  the  following  equations  for  ,  a^y ,  and  cr^  : 
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As  one  can  see  there  is  no  contribution  of  new  M-moments  in  the  equation  for  a^y .  The  equations  for  qy,  q^, 
ayy ,  a^^ ,  ay^  can  be  obtained  by  a  proper  rotations  of  indexes  in  Eqs.  (91)  -  (93).  Thus,  to  obtain  equations  for 
the  heat  flux  (qx,qy,qz)  and  the  stress  tensor  (cr^,cr^^,cr^2’‘^jy’‘^yz)  we  have  to  derive  equations  for  M-moments. 
Let  us  obtain  equations  for  Mi^iyiz  •  Multiplying  the  Boltzmann  Equation,  Eq.  (1),  by 
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we  obtain 


>-OT 
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(95) 


Transferring  the  terms  in  front  of  the  derivatives  inside  of  the  derivatives  brackets  in  Eq.  (95),  we  obtain 
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Substituting  Eqs.  (24)  and  (84)  into  Eq.  (96)  and  integrating  the  obtained  equation  over  the  entire  velocity  domain, 
we  obtain 
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Following  the  Chapman-Enskog  recipe  [5,  6],  let  us  put  r  ^  0  in  Eq.  (97);  then  substituting  Eq.  (58)  for  f^rad 
the  left-hand  side  of  the  obtained  equation  we  obtain  the  following  equation  for  (by  integrating  via 

Mathematica): 


M- 


Ixlylz 


-P-Vt  -[Xx-Xy-Xz-fu  -fl-d^X  =[^-^yz)-^  +  {^-^xz)^  +  [^-^xy)-  —  + 


16 


/  \  duj,  ,  ^  dwy  , 


)■— 4' 

’  V 


dt 

^  du^  du  ^ 


dt 


- 1 - 

dy  dx 


+  • 


f  du^  du^  , 

— -  +  — ^  +  (98) 

dz  dx 


16 


dUy  du^ 
dz  dy 


+  —  [4-Vf-ay,y-^i^4-Vf-a^^y-^i^4-Vf-ay^l. 


A 

dx ' 


Substituting  the  expression  for  dui/dt ,  Eq.  (15),  into  Eq.  (98)  we  obtain  the  final  expression  for  the  M 
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we  can  drop  terms  that  are  of  the  order  of  ;  in  the  right-hand  side  Eq.  (99).  This  yields 
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All  other  M-moments,  Eqs.  (86)  -  (89),  are  derived  in  Appendix  B.  After  collecting  all  M-moments  in  Eqs.  (91) 
and  (93)  together, 
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and  substituting  them  into  Eqs.  (91)  and  (93)  we  obtain  the  final  equations  for  and  ■  By  changing  the  order 
of  indexes  (x,y,z)  to  (y,x,z)  in  Eqs.  (91),  (102)  -  (104)  we  obtain  the  equation  for  qy ,  and,  similarly,  by  using 
(z,y,x)  the  equation  for  qy  ;  by  using  (x,z,y)  and  (y,z,x)  in  Eqs.  (92)  and  (105)  we  obtain  the  equations  for 
and  (jy^ ,  respectively;  and  by  using  (y,x,z)  in  Eq.  (92)  fwe  obtain  the  equation  for  <jyy .  Thus,  the  resulting  set  of 
13  equations,  Eqs.  (14)  -  (16),  (91)  -  (93),  (102)  -  (105)  with  rotations  of  indexes  as  described  above,  represents  a 
complete  set  of  the  13  regularized  Grad’s  moment  equations  for  p,  (ux,Uy,Uz),  Vj ,  ■> 

’  ^xy  ’  ^xz  ’  ^yy  ’  ^yz  )  • 

Finally,  the  coefficients  A14  -  A29  in  the  velocity  distribution  function  ,  Eq.  (31),  can  be  obtained  from 
Eq.  (86)  -  (90)  by  substituting  for  /  ;  which  yields: 


Ai4 


3S4-p-vf 


(106) 


Ai7  - 


^2x2y 

64-p-Vf  ’ 


^18  - 


^2x2z 

64-p-Vt  ’ 


Ai9  - 


^2y2z 

64-p-Vt  ’ 


(107) 


^20  -  ■ 


M 


3xly 


96-p-Vt 


^21  - 


3xlz 


96-p-Vt 


^22  -  ■ 


M 


3ylx 


96-p-Vt 


^23  -  ■ 


M 


3ylz 


96-p-Vt 


A24  - 


3zlx 


96-p-Vt 


^25  - 


M 


3zly 


96-p-Vt 


,  (108) 


A 


26  - 


^\x\y2z 

31-p-Vj  ’ 


A27  - 


^\y\z2x 

31-p-Vj  ’ 


^28  - 


^\z\y2x 

32-p-V}  ’ 


(109) 


A29  - 


^Ixlylz 

%-p-V^  ' 


(110) 


V.  CONCEUSIONS 
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We  have  presented  a  new  set  of  moment  equations  for  rarefied  gas  dynamics.  Our  equations  are  a  closure 
for  Grad’s  13  moment  equations  extended  to  the  third  order  of  the  Knudsen  number.  We  have  assume  a  Hermite 
polynomial  approximation  for  the  monatomic  gas  velocity  distribution  function,  the  BGK  approximation  of  the 
collision  term  in  the  Boltzmann  kinetic  equation  and  used  the  the  well-known  Chapman-Enskog  regularization 
method  that  has  been  previously  used  to  derive  a  closure  of  the  Euler’s  gas  dynamic  equations.  We  have  shown 
that  the  selected  29-term  Hermite  polynomial  representation  of  the  velocity  distribution  function  makes  good 
physical  sense  and  obtained  the  coefficients  for  this  polynomial.  Our  regularized  Grad’s  13  moment  equations 
differ  from  a  similar  set  of  equations  obtained  by  Struchtrup  and  Horrilhon  [4],  who  have  used  a  very  different 
and  complicated  method  and  have  not  assumed  a  Hermite  polynomial  representation  of  the  velocity  distribution 
function.  On  the  contrary,  the  closure  method  presented  in  this  paper  turns  out  to  be  quite  straightforward  and 
comprehensible. 

APPENDIX  A 

Let  us  obtain  the  general  moment  equation  for  .  Multiplying  the  Boltzmann  Equation,  Eq.  (1),  by 

+{Vy-Uyf  ,  (Al) 

we  obtain 

I  ■  -  M, )  ■  [(F,  -  a,)2  +  (F,  -  u,  f  +  (F,  -  a,  f  ]  ■  + 

+ 1  ■  -a,- )  ■  Fi  ■  [(F,  -  a,)2  +  (F,  -  a,,  f  +  (V,  -  a,  =  (A2) 

«yf  +  (K  f  ] ■  *(« ■  /)  ■ 

Transferring  the  terms  in  front  of  the  derivatives  inside  the  derivative  brackets  in  Eq.  (A2),  we  obtain 
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—  ^2  ^2^.  /■l_Lr .r/T/  ^2  ^2  ^2^.  /-Y 


^l^j-(^f-t^f)w,-u,r+(^y-uyr+(K-t^^n-fj+l^jW,-u,r+(^y-uyr+(K-t^^n-fj~+ 

+  “*  +{^V  +(^2  ““;■)■  (^A  ““A)'/’(^y  + 

+  M^-(p-(f^.-M,)-/-(^;-W;))-^  =  Y-(^-W/){(^x-^<x)^+(^y-^^j)^+(f^z-«z)^]-^^(«-/)- 


(A3) 


Integrating  Eq.  (A3)  over  the  entire  velocity  domain,  we  obtain 


A  A 

!(Vi  -«i)-[(v,-uyf  +(y,-«y^  +(v^-u,f]-/-<i^v 

y2.  y 


f 


X ■  J - "xf  +  (K. -Uyf  +  (K-u,f]-f-d^V 

V 


^  du  ^ 


dt 


P-\{Vi-Ui)-{Vj^-uj^)-f-d  V 


V  V 


dxt 


P 


y2  y 


•  [(Vt -Ui)-[{V, -u,f  +(E, -Uyf  +(V, -u,f]-{Vk -Uk)-f-d^V 


duk 

dt 


[(Vk -Uk){{V, -U,f  +(Vy -Uyf  +(V, -u,f]-f-d^V 


dxk 

'p 


\2  y 

^  P 
^k-- 
V  ^  V 

( 


u,f  +  (K,  - «,, )^  +  (('z  - "z )' ]  ■/  ■  << V 


“A-, 

V  ^  V 


SUj 

dxt 


•  l[(Vx -u,f  +(Vy -Uyf  +(V, -u,f]-f-d^V 


dxt 


p-\(Vi-Ui)-(Vk-Uk)-(Vj-Uj)-f-d^V 
V  V 

P_ 

2 


V 


J  J' 

p 


dj^ 

dxix 


Uk -Pf  {Vi-Ui)-{Vj -Uj)-f-d^V 
V  F  J 


dxK 


[(V,-u,)-  (V,-u,f+(Vy-Uyf+(V^-u^f  -Stin-d^V 


(A4) 


Expressing  the  third  and  ninth  terms  in  the  left-hand  side  of  Eq.  (A4)  as 


\  p-{Vi-Ui)-{Vk-Uk)-f-d^V 

\v 


d^ 

dt 


\P' 


V  V 


{Vi-Ui)-{Vk-Uk)-Sik-^ 


'N  ^ 

■f-d^V 

duty 
- -  + 

r 

P- 

2  J 

dt 

V 

2 

y 

QUj 

dt 
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^  ^  ■ 

Uk-p-\{Vi-Ui)-{Vj-Uj)-f-d^V 

V  y 


Uk-p-\ 

V 


y2  y 

(Vi-Ui)-(Vj-Uj)-dij-^ 


f-d^V 


Sxt 


Uk-P 


SUj 
dxk  ’ 


and  taking  into  account  Eqs.  (6)  and  (7)  and  the  fact  that  (jy  =  (jy ,  we  obtain 


dt 


5  p-Vf 


2^ 


2  2 

\  J 


dU; 


duu 


dU; 


-^  +  (^ik  +  ^ — \^k  -dij  +  dk  'T —  + 

dt  dt  dxj^  dxj^ 


5  p-Vp 

Uk - 

2  2 

J 


dU;  /  \ 


p-\iVi-Ui)-iVk-Uk)-iVj-Uj)-f-d^V 
V  K 


du  ; 

IT  =  (A5) 

dxh 


=  ^-\{Vi-Ui)-{V-uf-St{f)-d^V, 
2  ^ 


where 


t^2 

^ZZ  =  P-\  f  <V^-U^f  -di^V - ^=-CJyy-CJ^^. 

V  ^ 


(A6) 


Here  we  have  used  the  general  expression  for  stress  tensor  ay  given  by  Eq.  (7).  Substituting  duk  I  dt  from  Eq. 
(15)  into  Eq.  (A6)  and  after  some  simple  algebra,  the  equation  for  the  heat  flux  turns  into 


dq;  d  /  N  du; 

77  +  7 — -dij^dk’ - 
dt  dxj^  dxj^ 


^5  Vj^ 
J  ^  J 


d 


=  ^-\{Vi-u{)\V-uY-St{f). 

2  y 


-u^).(V-uf-(Vk -Uk)-f-d^V  +  p-\{Vi -UiYiVk -Uk)-(V,- -uP-f-d^V 


J 

\  f 


Sx,.  ^ 


P 


d 


<  p-Yp 


V 


da 


Jk 


dXj  J 


y2  y 


V  K 


J  J' 


dUj 

dxk 


(A7) 


.-r^2 


Next,  let  us  derive  a  general  equation  for  a^  .  Multiplying  the  Boltzmann  Equation,  Eq.  (1)  by 


fn-(V^-u^)-(Vy-Uy) 


(A8) 


we  obtain 


m-{V^-Uy)-iV^-Uy)- +m-(V^-Uy)-(V^-Uy)-Vk-  =m-{Vx-Uy)-{V^-Uy)- St(n  ■  /)  .  (A9) 


dxp 
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Transferring  the  terms  in  front  of  the  derivatives  inside  the  derivative  brackets  in  Eq.  (A9),  we  obtain 


0 

^[p-(V^-Ux)-{Vy-Uy)-f)+[p-{Vy-Uy)-f)--^  +  {p-{Vx-Ux)-f) 


dt 


dUy 

dt 


+  ^(p-(^k  -Uk)-{Vx-Ux)-{Vy-Uy)-  +  ■  P-(Vx-Ux)-{Vy-Uy)-f)  + 

OXk  OXj^ 

+  (p-{Vk-Uk)-{Vy-Uy)-f)-^  +  [p-Uk  ■{Vy-Uy)-f)-^  + 

OXk  ^^k 

+  {p-(Vk-Uk)-(Vx-Ux)-f)-^  +  {uk  ■p-(Vx-Ux)-f)-^  =  p-iVx-Ux)-iVy-Uy)-St(f)  . 

OXk  ^^k 


(AlO) 


Integrating  Eq.  (AlO)  over  the  entire  velocity  domain,  we  obtain 


P-\(Vx-t^x)-iVy-Uy)-f-d^V  +  p-\{Vx-Uy)-f-d^V 


V  V 
d 

H - 

dxk 

( 


V  V 


y  -yj 


dUx 

dt 


f 

P-  !(^x-^x)-f-d^^ 

V  V 


p-\iVk-Uk)-{Vx-Ux)-iVy-Uy)-f-d^V 
V  V 


dxi 


dUy 

dt 

Uk-p-\iVx-Ux)-{Vy-Uy)-f-d^V 

\  V  J 


p-\{Vk-Uk)-{Vy-Uy)-f-d^V 

y  V  j 

(  \ 

p-\{Vk-Uk)-{Vx-Ux)-f-d^V 

V  V 


dx^ 

dUy 

dxi. 


Uk-P-\iVy-Uy)-f-d  V 
V  J 


du^ 


Uk-p-\iVx-Ux)-f-(^^V 


X 

\  V  J 


dxk 

^  ^Uy 

dxk 


=  p-!(Fx-Ux)-(Fy-Uy)-St(f)-d¥  . 


(All) 


Expressing  the  sixth  and  eighth  terms  in  the  left-hand  side  of  Eq.  (A1 1)  as 


\  p-{Fk-Uk)-{Fi-Ui)f-d^F 

\v 


^  5m,- 

f  / 

_ l_  _ 

J  p- 

Fp 

■f-d^F 

duj 

- -  + 

f 

P- 

Fp 

2 

/ 

) 

dxk 

y 

2 

y 

du  i 

^,(A12) 

dX; 


where  indexes  ij  =  {xy,yx) ,  and  taking  into  account  Eq.  (4)  and  the  general  expression  for  the  stress  tensor,  Eq.  (7) 
we  obtain 


d  /  \  dUy.  dUy  Q 


t  \  ^^X 

~[‘^xy  )  + 


dt 


p-Fj 


dy  dx 


dxi 


(A13) 


p- i(Vk-Uk)-(Vx-Ux)-(Fy-Uy)-f-d  F  \  ^  p-  \(Fx-Ux)-(Fy-Uy)-St(f)-d  F  . 

y  V  J  V 


Now,  let  us  derive  a  general  equation  for  cr^  .  Multiplying  the  Boltzmann  Equation,  Eq.  (1)  by 


22 


(A14) 


we  obtain 


m  • 


2^ 


d{n-f) 

dt 


+  m' 


2^ 


Sxt 


^  7  Vt'^ 

V  y 


St(n-/).  (A15) 


Transferring  the  terms  in  front  of  the  derivatives  inside  the  derivative  brackets  in  Eq.  (A  15),  we  obtain 


r  f 


P- 

V  V 


2^ 


dxt 


f 

J  J 
f 


V  V  X  XJ  J  )  [  2  dt 


p-^k-Uk)- 


2  Vt 


2^ 


(v,-u,r- 


•/ 


dxt 


Hk -P- 


V 


2  kj" 


2^ 


iv,-u,r- 


f 


dUy 


+  (2.p.(Kj-«i).(F,-«,)./).^  +  (p.Ut.(K, -«,)./) 


5xi 

P<Vk-^k)  A  QVt  ,  i  H-P  ^ 


dxi 


dxK 


Dxk 


Stif)  . 


(A16) 


Integrating  Eq.  (A  16)  over  the  entire  velocity  domain,  we  obtain 


f  f 

p-\ 

\  v\ 


2^ 


f-d^V 


2-p-\{V,-u,)-f-d^V 


dUx 

dt 


P 


y2  y 


dVy 

dt 


dxk 

( 


P-\(Vk-Uk)' 

V 


2  Vy 


2\ 


iV.-Ux)- 


f-d^V 


dxt 


Uk-P-I 

V 


2  Vy 


2\ 


iv.-UxY- 


f-d^v 


2-p-\{Vy-Uy)-{Vx-Ux)-f-d^V 
V  V 


dxi 


P-^k- \(Vx-Ux)-f-d  V 
V  J 


dxt 


P 


\2  y 


■\iVk-Uy)-f-d^V 


dVy 

dxk 


Uk-p 


,3,7! 


^  2  y 


■ff-d^F  ■^  =  p-! 

dxk  V 


2  Vy 


2^ 


iVx-UxY- 


■Stif)-d^V. 


(A17) 


Using  Eq.  (A12)  in  the  sixth  term  in  the  left-hand  side  of  Eq.  (A17),  substituting  dVy  I  dt  from  Eq.  (16),  and  then 
taking  into  account  Eqs.  (3)  and  (4)  and  the  general  expression  for  the  stress  tensor,  Eq.  (7),  we  obtain 
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1  Tr2  duj  2  doj  2  du,-  d  (  \  ^  ^.2  ^^y 

£L.  p.r}.^.  (u,  cr„}+2^cr^t-^  +  p-rT'^  + 


dx;  3  dxi  3 


dxj  dx]^ 


+  ^\p-\\{Vy-UyY—^\-{Vk-Uk)-f-d^V\  =  p-\\{Vy-UyY-^\-St{f)-d^V. 


,2 


(A18) 


APPENDIX  B 


Let  us  obtain  an  equation  for  M2x2v  Multiplying  the  Boltzmann  Equation,  Eq.  (1),  with  the  BGK 


approximation  substituted  for  the  collision  term,  Eq.  (24),  by  m  ■  (4  ■ 


(Fx-UxY-2 


■vi\{^-(Vy-UyY  -2-vi)  we 


obtain 


+  m-Vi,-(4-(V,-up -2-v}y(4-(V,-up -2-r})-?YY  = 

^m-(4-(V,-up -2-v})-[4-(y,-up -2-vf  )-n-dM^  . 


Transferring  the  terms  in  front  of  the  derivatives  inside  the  derivatives  brackets  in  Eq.  (Bl),  we  obtain 


j^{^4-p-(2-(r,-u,r-r})y2-(r^-u,f-v})-fy[p-s-((v,-u,f  +  (v,-u,f-rf)-fy 
+is.p.(v^-u,)y2-(r,.-up-v}yfy^+is-p-(ry-u,)\2-(v,-up-v}yfy?Y+ 


+-^[4-piVt-uty(2-{V,-up-Vfy[2-(V^-up-Vr^yfj+ 

+  ^(4'P'%'(2.(K,-u,)2-Kf)-(2.(K,-a,.)2-Fr4/)  + 

-  -  ~  ~ 

+%■  p-{Vj^-Uk)-  (^(Fx-UxY+iFy-UyY-Frj-f  ((^X  ~^x)^  +(^y  ~^y)^  ~^Tj'f  ■ 

+16-p-(Vk-Uk)-j^(Vx-Ux)-(2-(ry-UyY-Vr^yj^  +  16-p-Uk-j^(Vx-Ux)-(2-(ry-UyY-rr^y  ■ 

+16-p-(Vk-Uk)-l(Vy-Uy)-l2-(rx-UxY-rr^)fj~  +  i6-p-Uk-l(Vy-Uy)-l2-(rx-UxY-rr^)f  • 
=  ]..p.j^4-(rx-UxY-2-Vr^j-(4-(ry-UyY-2-Vr^j-(/j^-/)  . 
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Following  the  closure  method  described  in  Section  IV  for  calculating  moment  ,  we  substitute  f^rad  / 

into  the  left-hand  side  of  Eq.  (B3)  and  fcrad  +  ^’  fu  '  f\  i^^to  the  right-hand  side  and  then  integrate  the  resulting 
equation  over  the  entire  velocity  domain  using  Mathematica.  This  yields 


Substituting  duf  /dt  and  dV^  Idt  from  the  momentum  and  energy  conservation  law  equations,  Eqs.  (15)  and  (16), 
into  Eq.  (B4)  we  obtain 


32  q,-VT^ 
5  p 


dp 

dx 


^32  qy-Vf^ 


dp 

dy 


+  1  —-Vj 

5 


dqy.  ( 32  2 

- Vt 

dx  15 


64 


dV^  f  64 


16 


+  y.ri?.(2. 


J 

dv}  f32  )  dV^  16  ,0  L 


Sqy 


dUx 

dx 


(j^y.  -  a 


\  (JUy  16  ( 

r  dy  3  -{^xx  +  ^^yy, 


32-V^ -a 


xy) 


\(dUx  ^  S^y 

'  dy  dx 


32  o-,,  32  o- 


V 


■  +  ■ 
p  3 


P 


32  o-,,  32  o- 


■  +  ■ 
p  3 


P 


Sqk 

dxk 

'^ik 


5  p 


da, 


kx 


dxu 


5  p 


da 


ky 


dxu 


duj 

dxy 


M 


2x2y 


(B5) 


Taking  into  account  Eq.  (100)  and  dropping  all  terms  of  order  in  Eq.  (B5)  we  obtain. 


^2x2y  32  Vf  (  dp  dp 

T  5  /9  i  dx  ^  dy 


+  1  --V^ 
5 


32 

+  ly^z 

+  ^-F^(2-, 


dVj  ( 64 


dV^  ( 64 

— - q 

dx  I  5 


XX  yy } 


dUy 

y 


^qx  1 

f  32 

•Erl 

dqy 

dx 

b ' 

) 

1  dV^ 

16 

+  — 

■V^  • 

(2  •  a,,,,  -  c 

'xx)‘ 

du^ 

— -  + 

3 

V  yy 

dx 

1  * 

du^ 

+  (32. 

•  Vy  •  Cyy  )• 

''  dUx 

dUy 

yy  / 

dz 

\ 

1  xy  J  \ 

V  ^ 

dx 

(B6) 


Changing  the  order  of  indexes  (x,y,z)  to  (x,z,y)  in  Eq.  (B6),  we  obtain  an  equation  the  M^xjz  moment,  and, 
similarly,  for  the  M2y2z  moment. 
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Let  us  obtain  an  equation  for  .  Multiplying  the  Boltzmann  Equation,  Eq.  (1),  with  the  BGK 
approximation  of  the  collision  term,  Eq.  (24)  by  ~Ux)^  -48-F^  -(V^  ~Ux)^  +12-F^]  we  obtain 


4\ 


■  i\6-(y,-u,f  -4S-yf  -(v.-u.f  +n-vf) 


4\  d{f-n) 


dxy 


(B7) 


=  m 


-4S-v^  ■(V.-U.f  +n-yf)-„ 


fu  ~f 


Transferring  the  terms  in  front  of  the  derivatives  inside  the  derivatives  brackets  in  Eq.  (B7),  we  obtain 


P 


+{^p\(>4-(y,-up-<>6-y}4v,-u,)\f]-^+-^(p4yi-ut)\\(>4y,-up-4i-y}4v,-u,f+n-v^\fy 
+-^[p-H\i(,-{v,-u^f  -4ii-y}  -(y.-up  +  n-v^\ f)  +  p-iy,-ut)-[[4t,4v,-u,f  -24-y}\ f 


dVj 
J  dxy 


(B8) 


+p-Uj^- 


+p-Uj^- 


(48-(F,-M,)2-24-FrT/ 


BVt 
J  dxy 


+p- 


(Vt-Uk){(A4y,-up -(V.-u,))-/ 


[64-(v,-up-96-yP{v,-u,))-fj-^^p-(i6-(y,-up-4!i-yP{y,-up+u-y^)- 


dUy. 

I 

J  dxj^ 

4  \  fu  ~f 


Following  the  same  procedure  as  for  the  moment  M2x2y  described  in  the  paragraph  following  after  Eq.  (B2)  we 
obtain 


M, 


4x 


384 


du 


dVf  (480 


= - qx - ^  +  (48-cr^) - h - ^^+48-m^-c7- 


dt 


dt 


Yitj 

''^P\ 

)  dx 

5 

dqx 

dx 


ttI 

384  ) 

dUy  { 

'384  ) 

dUx 

(384  ) 

dUx 

^xx'^T  + 

— -“x'^x  |- 

dx  i, 

-■Uy-^lxV 

dy 

+  1  •  “z  •  ^x  1 

dz 

(  96 

+  48-m^  -cr^ 

)  dV^ 

r96 

(  dv^ 

J'lT 

+ 1  —  •  +  48 

•Mz 

(B9) 


Substituting  duj  I  dt  and  dVj  Idt  from  the  momentum  and  energy  conservation  law  equations,  Eqs.  (15)  and  (16), 
into  Eq.  (B9)  and  then  dropping  all  terms  of  the  order  in  the  resulting  equation,  we  obtain 
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(BIO) 


Max 

T 


Sqx 

dx 


192  ■qx-Vj  '\  dp 

5  p  dx 
V  ^  J 


(32 


V 


du^  SUy 
dx  dy 


du^'] 


Changing  the  order  of  indexes  (x,y,z)  to  (y,x,z)  in  Eq.  (BIO),  we  obtain  an  equation  the  M^y  moment,  similarly, 
for  the  Maz  moment. 

Let  us  obtain  an  equation  for  M^xiy  •  Multiplying  the  Boltzmann  Equation,  Eq.  (1)  with  the  BGK 
approximation  of  the  collision  term,  Eq.  (24)  by  m-8-^2-(F^  -1-Vj  -(1^  ~^y)  obtain 


S-n,-(2-(V,-uy  + 

+!i-m-ri-[2Ar,-uy-i-vP(v,-u,))-{v,-Uy)-q^^  (Bll) 

=  i-m-[2-(V,-uy-3-Vy{V,-iO)iVy-Uy)-n-qM-q-  . 

Transferring  the  terms  in  front  of  the  derivatives  inside  the  derivatives  brackets  in  Eq.  (Bll),  we  obtain 

2 

+(8.p.(2.(F,-u,)3-3.Kf  (F,-a,))./).^  +  (24.p.(2.(K,-uy-»f).(^',.-»,.)-/)'^+ 

+  £-{»-P<Vi-Ui)-(2AVy-uj-3-VyiVy-Uy)){Vy-Uy)-f  y 
+^(8'/2-ai-(2.(F,-a,)3-3.Fi?.(F,-u,)).(F,-Uj,)./)+ 

(B12) 

+(^g- p-(Vt  -u^)-(2-(V,-up -l-vPiVy-uy)- f  +  p-ut  ■(2-(V,-uy -3-vy(Vy-uy'j-f  'j-^  + 

+  {^24-p-(Vt-Ut){2-(Vy-Uy-V.^-y{Vy-Uy}f)-^+[24-p-Ut-{2Ary-Uy-V;-}{Vy-Uy}f)-^  = 

=  %■  p-{2AVy-Up  -3-Vy(Vy-Uy)\(Vy-U,yq^i^  . 

Following  the  same  procedure  as  for  the  moment  M2x2y  described  in  the  paragraph  after  Eq.  (B2)  we  obtain 


^-rrl. 

+2A-P-  (Vk-uj^)-(Vx-Ux)\Vy-Uy)-f  p-uj^  -  (Vx-Ux)-(v^ 


y  Uy 


dVy 

dxp 
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Substituting  duj  /dt  and  dV^  Idt  from  the  momentum  and  energy  conservation  law  equations,  Eqs.  (15)  and  (16), 
into  Eq.  (B13)  and  then  dropping  all  terms  of  order  in  the  resulting  equation  we  obtain 
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(B14) 


Changing  the  order  of  indexes  (x,y,z)  to  (x,z,y)  in  Eq.  (B14),  we  obtain  an  equation  the  M2xiz  moment;  for 
(y,x,z)  we  obtain  an  equation  the  moment,  and  similarly  for  the  »  ^3ziy  ’  M32i;^  moments. 

Finally,  let  us  obtain  an  equation  for  •  Multiplying  the  Boltzmann  Equation,  Eq.  (1)  with  the  BGK 

approximation  of  the  collision  term,  Eq.  (24)  by  m-8-(V^  -u^)^  we  obtain 


+  8-m-Vf(Fy-Uy)-(Vy-i,y)-(2-(K-uy-V^')-yP^^  (B15) 

=  8-m-m-8-(Fy-u,)-yy-Uyj-(^2-(Vy-Uyf-Vy'j-rrP*^  ^  . 
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dVj 

dt 
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(B16) 


Following  the  same  procedure  as  for  the  moment  M2x2y  described  in  the  paragraph  after  Eq.  (B2)  we  obtain 
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(B17) 


-(l6-cT^z-Fr)-^  +  (^y-^x-“^]-^  +  (l6-cT^z-Fr)-^  +  (l6-cT^z-Fr)-^  +  (l6-cT^-Fr2). 


du^ 

dz 


Ml 


Ixlylz 


Substituting  du, !  dt  and  dVj  Idt  from  the  momentum  and  energy  conservation  law  equations,  Eqs.  (15)  and  (16), 
into  Eq.  (B17)  and  then  dropping  all  terms  of  the  order  in  the  resulting  equation  we  obtain 
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(B18) 
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Changing  the  order  of  indexes  (x,y,z)  to  (x,z,y)  in  Eq.  (B18)  we  obtain  an  equation  for  the  moment, 

and  using  (z,y,x)  we  obtain  an  equation  for  the  moment. 
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